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Accurate and efficient online parameter identification and state estimation are crucial for leveraging Digital 
Twin simulations to optimize the operation of near-carbon-free nuclear energy systems. In previous studies, we 
developed a reactor operation digital twin (RODT). However, non-differentiabilities and discontinuities arise 
when employing machine-learning-based surrogate forward models, challenging traditional gradient-based in- 
verse methods and their variants. This study investigated deterministic and metaheuristic algorithms and devel- 
oped hybrid algorithms to address these issues. An efficient modular RODT software framework that incorpo- 
rates these methods into its post-evaluation module is presented for comprehensive comparison. The methods 
were rigorously assessed based on convergence profiles, stability with respect to noise, and computational per- 
formance. The numerical results show that the hybrid KNNLHS algorithm excels in real-time online applica- 
tions, balancing accuracy and efficiency with a prediction error rate of only 1% and processing times of less 
than 0.1 s. Contrastingly, algorithms such as FSA, DE, and ADE, although slightly slower (approximately 1 s), 
demonstrated higher accuracy with a 0.3% relative L2 error, which advances КОРТ methodologies to harness 
machine learning and system modeling for improved reactor monitoring, systematic diagnosis of off-normal 
events, and lifetime management strategies. The developed modular software and novel optimization methods 
presented offer pathways to realize the full potential of RODT for transforming energy engineering practices. 
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I. INTRODUCTION 
A. Concept of Digital Twins (DTs) 


The concept of digital twins (DTs) has gained signifi- 
cant traction in recent years, particularly in the engineer- 
ing and industrial disciplines [1]. A DT refers to a virtual 
model that closely mirrors a real-world physical prod- 
uct, system, or process, serving as an effectively indistin- 
guishable digital counterpart for practical purposes such 
as simulation, integration, testing, monitoring, and main- 
tenance. DTs have been recognized as the fundamental 
premise of product lifecycle management, encompassing 
the entire lifecycle of the physical entity they represent, 
including creation, construction, operation/support, and 
disposal [2—6]. 

The value and fidelity of a DT representation depend 
on the specific use cases, considering that its granularity 
and complexity are tailored accordingly. Notably, a DT 
can be conceptualized and utilized even before its physi- 
cal counterpart exists, allowing for comprehensive mod- 
eling and simulation of the intended entity's lifecycle. 
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The origin of the DT concept can be traced back to 
NASA, which pioneered its practical definition in 2010 
as part of its efforts to enhance the simulation of space- 
craft using physical model representations [2, 6]. Since 
then, the development of DTs has progressed in tan- 
dem with advances in product design and engineering. 
Traditional manual drafting techniques have given way 
to computer-aided drafting and design and model-based 
system engineering approaches, which establish a direct 
link between the DT and signals from its physical coun- 
terparts. 

In recent years, DTs have garnered significant atten- 
tion in various industrial sectors. Notably, their appli- 
cation in the chemical processing industry has demon- 
strated versatility [7]. In addition, DTs are increasingly 
utilized in simulation-based vehicle certification and fleet 
management [8], as well as in addressing the challenges 
of unpredictable and undesirable behaviors in complex 
systems [9, 10]. A specific area of growth has been ob- 
served in the nuclear plant industry, where digital twin 
architectures have been developed for the management 
and monitoring of nuclear plants, reflecting the expand- 
ing scope and potential of technology [11, 12]. 


B. Previous Work on Reactor Operation Digital Twins 
(RODTs) 


A significant amount of research has been conducted 
to explore the potential applications of reactor opera- 
tion digital twins (RODT) in the nuclear energy domain. 


КОРТ represent a specific instantiation of the DT con- 
cept, focusing on the numerical representation of nuclear 
reactors for real-time prediction, optimization, monitor- 
ing, control, and decision-making during their opera- 
tional stage [13]. 

A previous study [13] introduced the first prototype of 
an RODT specifically tailored for the prediction of neu- 
tron flux and power fields in the operational stage of the 
HPR1000 reactor core [14]. The prototype demonstrated 
the feasibility of using the RODT for online monitor- 
ing and prediction of key reactor parameters, thereby en- 
hancing operational understanding and performance as- 
sessment. 

In [15], the forward model of the RODT was realized 
through the application of a nonintrusive reduced-order 
model constructed using an SVD autoencoder to learn 
the nonlinearity of the field distribution. Machine learn- 
ing techniques, specifically leveraging k-nearest neigh- 
bour (KNN) and decision trees, are employed to estab- 
lish forward mapping. Additionally, the inverse model 
adopts a generalized latent assimilation method for ac- 
curate estimation of the model parameters. 

While demonstrating promise, opportunities remain 
for improving the efficacy of RODT for practical deploy- 
ment. In particular, the inverse solver relies on a compu- 
tationally expensive methodology that limits its potential 
integration within an online framework. Therefore, the 
aim of this work [16] is to enhance the key capabilities of 
the RODT. This includes proposing an advanced differ- 
ential evolution algorithm to upgrade the inverse solver 
for improved efficiency and accuracy of parameter esti- 
mation. In addition, uncertainty quantification was con- 
ducted to validate the RODT's performance considering 
noisy observational data. Numerical validation experi- 
ments were performed across representative domains for 
an HPR1000 pressurized water reactor core to demon- 
strate its potential for engineering applications. 

The key steps for constructing an RODT include: (1) 
training a reduced-order model (ROM) using model or- 
der reduction (MOR) [17—20] techniques as well as an 
efficient forward model and (ii) adaptation and imple- 
mentation of inverse models to infer the input parameters 
and the resulting field distribution. 


C. Introduction to Reduced Order Models (ROMs) 


Many mathematical models encountered in real-life 
processes present computational challenges when em- 
ployed in numerical simulations because of their inherent 
complexity and large dimensionality. To address these 
challenges, MOR techniques have been developed to re- 
duce the computational complexity of these problems, 
particularly in the context of simulations involving large- 
scale dynamic and control systems. The ROM was com- 
puted as an approximation of the original model by re- 


ducing the dimensionality or degrees of freedom associ- 
ated with the model. 

The preparation of an ROM model involves two dis- 
tinct phases, the offline and online phases, which are 
integral to the creation of a simplified representation 
of a high-fidelity full-order numerical model, enabling 
efficient real-time simulations and analysis. The of- 
fline phase focuses on determining the inherent low- 
dimensional structure of the underlying full-order model 
and deriving the reduced basis functions from high- 
fidelity snapshots of the system. Intrusive MOR tech- 
niques [21] are employed when a comprehensive un- 
derstanding of the governing equations and numerical 
strategies employed in the full-order model is available. 
This requires access to the detailed numerical framework 
of the full-order model, which may involve proprietary 
commercial codes. By projecting a full-order model onto 
a reduced space, intrusive ROMs provide a concise rep- 
resentation with reduced computational complexity. 

However, in practical engineering scenarios [22—24], 
detailed information regarding a full-order model is of- 
ten unattainable, and the solvers or codes implement- 
ing these models are treated as black-box entities. This 
limitation poses a challenge to the implementation of 
traditional intrusive MOR methods. To address this is- 
sue, nonintrusive MOR techniques have been developed. 
Non-intrusive MOR approaches aim to construct ROMs 
without relying on detailed knowledge of the numerical 
framework of a full-order model. Instead, these meth- 
ods establish an input-output mapping between the in- 
put parameters and the reduced basis through data-driven 
techniques such as interpolation, regression, or machine 
learning, approximating the behavior of the full-order 
model based on a reduced set of parameters and basis 
functions. 

ROMs prepared through the offline and online phases 
play a critical role in the development of RODTs. 
RODTs require the integration of high-fidelity mod- 
els and solvers with real-time capabilities, and ROM- 
based emulators provide promising solutions [21—23, 
25]. Leveraging ROM-based forward solvers, RODTs 
can perform both one-to-one forward simulations and 
one-to-many simulations for inverse problems, encom- 
passing parameter identification [26—29], data assimi- 
lation [24, 30-35], sensitivity/uncertainty quantification 
[36—40], and optimization/control [37, 41-47]. 


D. Forward Model Construction 


In line with the framework RODT, as investigated in 
a previous study [13, 15], we selected KNN [48] to con- 
struct the forward model. The choice of KNN is justified 
by its numerous advantages, which render it a suitable 
candidate for the forward model in this context. 

First, it is a nonparametric algorithm, meaning that it 


does not assume any specific functional form for the rela- 
tionship between the inputs and outputs. Second, KNN is 
relatively simple to implement and computes efficiently, 
particularly for low-dimensional input problems. The 
complexity of the algorithm depends primarily on the 
number of training samples, making it suitable for prob- 
lems with moderate dataset sizes [49]. In a recent study 
[50], an investigation was conducted to explore the pos- 
sibility of utilizing different values of K for each dataset, 
considering the information provided by the correlation 
matrix. 


E. Inverse Model Construction 


Generally, the inverse problem can be divided into two 
categories: (1) the first category involves solving the for- 
ward system with knowledge of an existing dataset of in- 
put and output; (11) the second category involves solving 
the system's input with knowledge of the forward sys- 
tem and output [51]. In our case, the first category corre- 
sponds to the construction of the forward model, where 
we build a non-intrusive model as a surrogate to the full- 
order system. In the second category, we primarily focus 
on two specific scenarios: 


* Parameter Identification. Solving the parame- 
ter vector based on the information from observa- 
tion vectors: In this scenario, we aim to estimate 
the unknown parameters of the system based on 
the available observation vectors. The observation 
vectors contain information about the system's be- 
havior or response under conditions that are char- 
acterized by the parameter vector. 


State Estimation. Solving the field distribution 
based on the information from observation vec- 
tors: In this scenario, the objective is to determine 
the unknown field distribution of the nuclear sys- 
tem based on the available observation vectors. 


F. Contribution of this Work 


A previous work [13, 15, 16] demonstrated the feasi- 
bility of employing the RODT framework in engineering 
applications. However, a need exists for a systematic and 
in-depth investigation of the core components of RODT, 
namely parameter identification and state estimation. In 
addition, the development of mature and modular soft- 
ware is crucial to facilitating the application of RODT 
in practical scenarios. This motivated us to develop an 
efficient and appropriate surrogate inverse model to ac- 
curately predict the parameters and states. This also mo- 
tivated our research on hybrid optimization approaches 
that can effectively balance coarse- and fine-grid search 


(coarse-fine-grid search) methods, which are discussed 
in the following sections. 
The contributions of this research are as follows: 


° Investigation and implementation of novel 
gradient-free optimization algorithms specifically 
tailored to handle the challenges posed by the 
discontinuous surrogate forward mapping con- 
structed with the KNN forward model of our 
RODT and comparing the convergence profile, 
stability, and accuracy-time performance. 


Exploration and evaluation of the proposed 
coarse-fine-grid search methods within various 
global optimization approaches, aiming to strike 
a balance between accuracy and computational ef- 
ficiency in parameter identification and state esti- 
mation. 


* Development and deployment of efficient modu- 
lar RODT software integrated with the aforemen- 
tioned algorithms into the inverse problem solver, 
which applies to other sustainable energy systems. 


These contributions focus on advancing state-of-the- 
art RODT by addressing the unique challenges encoun- 
tered in solving inverse problems in nuclear engineer- 
ing and proposing novel methodologies for improved pa- 
rameter identification and state estimation. Furthermore, 
our work addresses the critical challenge of achieving 
real-time, high-precision simulations in nuclear energy 
DT applications, significantly enhancing the capabilities 
of DTs in the nuclear energy sector and facilitating im- 
proved monitoring, control, and optimization of nuclear 
energy systems. 


I. METHODOLOGY 


The modeling process of an RODT is classified into 
offline and online stages. 


* Offline phase: training a non-intrusive ROM, 
which includes (i) preparation of the full order 
snapshots; (ii) preparation of the reduced basis; 
(iii) preparation of the coefficient; and (iv) prepa- 
ration of the surrogate forward model to predict the 
coefficient with the information of the parameter. 


Online phase: parameter identification and state 
estimation, which includes (1) building an inverse 
model to predict the parameter based on the infor- 
mation of the observation data (clean or with ob- 
servation noise) and (ii) building an inverse model 
to predict the coefficient based on the information 
of the observation data (clean or with observation 
noise). 


The methodology is illustrated in the following 
flowchart. 
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Fig. 1: Working flow of the RODT schema. 


A. Offline phase: Training a Non-Intrusive ROM 


The variation in the power field ® in a nuclear core 
is characterized by physical laws that are generally de- 
scribed implicitly by a governing equation [52], such 
as the neutron transport or diffusion equations, which is 
written as Equation 2.1. 


F(®(r, р), p) = 0 
reOQcR^4 , 
нє? с КР 


Q.1) 


where Q C В“ represents the spatial domain of dimen- 
sion d, with d > 1, and Ф(и) = (r, p) is defined in a 
Hilbert space that is equipped with an inner product (-, -), 
and the induced norm ||.|| = ./(-,-). D C R? represents 
the p-dimensional feasible parameter domain that covers 
the operation of the reactor core. In this study, Eq. 2.1 
represents two-group neutron diffusion equations [53]. 


1. Preparation of the Full-Order Snapshots 


In this study, the numerical solution of Eq. 2.1 was 
obtained using the proprietary CORCA-3D code pack- 
age developed at the Nuclear Power Institute of China 
(NPIC) [54]. The CORCA-3D code is a black-box 
solver; it can be substituted with other proprietary codes 
as needed. 


To construct a non-intrusive ROM, it is necessary to 
gather a set of full-order solution snapshots. We sample 
ш: in the parameter domain D to obtain the discrete set 
D = {шщ € R? | pw; € 2,1 = 1,..., P), representing 
the entire set D. Using CORCA-3D, we obtain a solution 
snapshot set M = (ó(u) € RY | u € D). 

The knowledge of this manifold is implicitly depen- 
dent on the knowledge of the parameters jz, and for bet- 
ter construction of a non-intrusive model, we endeavor to 
find a low-rank representation for the solution snapshot 
set M. Reduced-basis methods (RB) [55-57] 

indicates that in state where Ф is sufficiently regular, it 
can be approximated by an n-dimensional reduced basis 
{qi}, as follows: 


Ф(и) € ó, (и) = 2. ou(u)qi- (2.2) 


For simplicity, we denote the approximation of the so- 
lution manifold as follows: 


ó»(u) = Qa(p) , (2.3) 


where Q € КУ х" represents the assembled basis N x n- 
dimensional matrix {q;, i = 1,...,n}, and o(u) € R” 


denotes the assembled coefficient n-dimensional vectors 
{aii = 1, eH, n}. 


2. Preparation of the Reduced Basis 


For simplicity and time efficiency, we chose the sin- 
gular value decomposition (SVD) method for generating 
the reduced basis. The steps are as follows: 

Correlation Matrix. The correlation matrix C is 
computed by taking the inner product of each pair of 
snapshots, as expressed by the following equation: 

Cig = 20008), A), VISEISP. Q4 

Here, (d(x), 6(p5)) represent the inner product of 
the ith and jth snapshots. Subsequently, the eigenvalues 
A; and corresponding eigenvectors v of the correlation 
matrix C are computed. 

Proper orthogonal decomposition (POD). The jth 
POD basis vector qJ, which is independent of parameter 
LL, is obtained as a linear combination of snapshots. 


N . 
q = wie. (2.5) 


i=l 


In this equation, v? represents the ith element of the 
j-th eigenvector, and ф* corresponds to the i-th snapshot. 
The magnitude of the j-th eigenvalue A; provides infor- 
mation about the relative importance of the jth POD ba- 
sis vector. 

Collection of the best basis minimizing the Frobe- 
nious error. The solution snapshot set is denoted as 
M = {¢(ux)}h_, C R^. When applied to computa- 
tion, we define the snapshot matrix as S Є RN*P. which 
contains snapshots ¢ (uk) as columns. The first n POD 
bases are assembled into a matrix Qn = [q1,-.--, qn] € 
Вх". Among all orthonormal bases of size т, ће POD 
basis minimizes the Frobenious norm least squares error 


m n 2 
V 22i 25j leül ) of 
the reconstruction of snapshot matrix S, which further 
corresponds to the Echart- Young theorem [58]: 


(defined in the form of || e |F = 


N 
mm|[S-Q.QrS|,- Уу ^ 06 


k=n+1 


where A; represents the kth singular value of matrix S. 
In summary, the POD basis constitutes a set of orthonor- 
mal vectors that offers the most effective n-dimensional 
representation of the given snapshots. 


3. Preparation of the Coefficient 


Following the aforementioned procedure, the coeffi- 
cient is constructed as the projection of the snapshot ma- 
trix S onto the vector space spanned by the POD basis 


Ѕрап{ д, q2;..- T 


For any (jt) belonging to the set S, the n- 
dimensional approximation of Ф( p) is given by Equation 
2.2. In this manner, the coefficient can be computed us- 
ing an orthogonal projection as follows: 


aq) = у $0295 - S qa), ол 


i=l (qi, qi 


where (q;, qi) equals 1 owing to the orthogonormality of 
the POD basis. 

Moreover, in matrix form, the coefficient matrix A, = 
[oa (и), ..., o. (u)]* € В"? follows the following re- 
lationship: 


S S Sn =QnAn=QnQzSs. (2.8) 


4. Training a Forward Model 


The goal of the forward model is to construct an ef- 
ficient surrogate mapping model between the parameter 
and coefficient, as follows: 


a(n) = ЕМ (и), 


where FME denotes the surrogate forward mapping 
constructed by the machine learning approach, such as 
model-dependent machine learning methods or model- 
free optimization searching methods. 

The response of the ROM model, denoted as FM" (y), 
corresponds to the predicted coefficient data. The image 
of space D under the forward map М (p) represents 
the set of responses for all possible states. The difference 
Y, — HQF™" (ш) is referred to as the observation data 
misfit or the residuals associated with parameter p. 

We design the cost function as Equation 2.12 based 
on the residual. The forward model is used in the 
search of the optimal approximation for the operator ¢ in 
the functional set (QFML FML є £(D, R")), where 
€ (D, R”) denotes the set of functions mapping from D 
to R”, with n the reduced dimension in ће MOR pro- 
cess. 

Previous research [13] evaluated various forward 
model approaches, among which the KNN algorithm 
demonstrated exceptional accuracy and efficiency. In 
this study, we specifically selected the KNN algorithm 
[48] as the forward model to exploit its full potential. 
To enhance the performance of the KNN algorithm, we 
focused on two crucial parameters that significantly in- 
fluence its efficacy. 

The first parameter relates to the voting distance in 
KNN, wherein we investigate the impact of utilizing ei- 
ther the Euclidean distance or Manhattan distance as 
the metric for determining the nearest parameter choice 
in the training set. Notably, these two distances can 
be expressed as Minkowski norms with an exponent of 


(2.9) 


p — 1,2, corresponding to the Euclidean and Manhattan 
distances, respectively. The Minkowski distance, a mea- 
sure of the discrepancy between the testing and training 
parameters, can be mathematically represented as fol- 
lows: 


dim(p) 
test train|p 
> | ИИ ; 


ї=1 


test __ кон __ 
p 


|u” — ш 

(2.10) 
where dim(s) denotes the dimension of the parameter 
ш. Finally, we choose the Euclidean metric. 

Another vital parameter is K, which refers to the num- 
ber of candidates voted for. The optimal choice of K can 
be determined by plotting the relative Lə error of the pre- 
dicted observation and true vectors, which is expressed 
as follows. 

From Figure 2(a), we can deduce that when K = 5, our 
forward model using KNN performs with the best predic- 
tion accuracy, and the variation in the error is relatively 
smooth, which corresponds to the continuous variation 
in the value of the field at the corresponding position of 
the field manifold snapshot. 

In contrast, the stair-shaped observation data predic- 
tions in Figure 3 demonstrate that the KNN forward sur- 
rogate model has a step-like progression as neighbors are 
adjusted between discrete datasets and thus are not con- 
tinuous or differentiable. In this manner, well-developed 
state-of-the-art algorithms based on gradient information 
for solving optimization and inverse problems did not 
work well in our case. To meet this demand, we inves- 
tigated other novel approaches, which are discussed in 
detail in the following sections. 


B. Online phase: Parameter Identification and State 
Estimation 


The inverse problem in nuclear engineering refers to 
the task of determining the model parameters that yield 
the observed data, as opposed to the forward problem, 
which involves predicting the data based on the given 
model parameters. The objective is to determine the 
model parameters 44* that approximately satisfy Equa- 
tion 2.12. 

The nature of the inverse problem depends on whether 
operator is linear or nonlinear. In most cases, particu- 
larly when solving systems governed by the neutron flux 
diffusion function, the inverse problem is nonlinear, ow- 
ing to the nonlinearity of the forward map. Thus, we at- 
tempt to develop nonlinear algorithms, gradient-free al- 
gorithms, or algorithms based on the initial guess given 
by fast linear algorithms, such as KNN. 
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Fig. 2: Finding the optimal K for the KNN model 
(Forward on the left, Inverse on the right). 


l. Setup of Inverse Model and Optimization 


In the context of optimal control theory, the governing 
equations that describe the behavior of a physical system 
are commonly referred to as state equations. However, 
in many practical scenarios, interest lies not only in the 
physical state itself but also in its impact on certain ob- 
jects or quantities. Furthermore, only a limited amount 
of data can often be obtained from the physical state. To 
address these considerations, an additional operator de- 
noted by H, known as the observation operator, is in- 
troduced. This operator maps the state of the physical 
system denoted by 6 to the desired observations denoted 
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Fig. 3: Non-differentiable stair-shaped observation data 
prediction via the KNN surrogate forward model (tested 
when the first parameter 5, varies from 0 to 615 steps, 
the other parameters take the midpoint values of the 
parameter intervals, respectively, 1250 for В„; 50 for 
Р„; and 295 for Tin). 


by Y, € R”, where m denotes the total number of ob- 
servations. Given a field $, the observation vector Y, is 
given by the following equation and approximated by an 
error: 


Y, = H? = HQF™" (u) + e. (2.11) 


where e € R” represents the observation noise vector 
that captures the presence of observational errors. The 
quality of a given set of model parameters pz is evaluated 
using the Lə distance between the simulated and practi- 
cal observations, which is quantified using the following 


cost function: 


2 
(и) := ||nQF ^ (p) – Yoll;, - Q.12) 

The objective of the inverse problem is to minimize 
the observation error, which can be formulated as 


p^ = argminJ(p) :— argmin (Inr qu) Е a ) : 
ncn | 


HED 
(2.13) 
It is important to note that the cost function described 
above is possibly nonlinear and non-differentiable with 
respect to the discontinuous function FM, Therefore, 
we developed the following algorithms to solve the opti- 
mization problem: 


HI. COMPARISON OF DIFFERENT APPROACHES 
FOR SOLVING INVERSE MODELS 


To solve inverse problems, machine learning and opti- 
mization methods can be employed to search for the op- 
timal parameters or field distributions that best fit the ob- 
served data [59]. Additionally, as mentioned in Section 
ILA 4, when we employ the KNN algorithm to construct 
the surrogate forward model, the surrogate forward map- 
ping is non-differentiable and discontinuous. This poses 
a challenge because the traditional approaches used to 
solve continuous optimization and inverse problems [60] 
are not applicable in this context. This necessitates the 
exploration of novel inverse problem-solving methods 
suitable for our particular problem. Hence, this study in- 
vestigates various deterministic optimization algorithms, 
metaheuristic algorithms, and their hybrids to address 
this challenge. 


A. Exhaustive Direct Search with Latin Hypercube 
Sampling (LHS) 


To evaluate and compare the different approaches to 
solving inverse models, it is essential to establish a 
benchmark for measuring their performance. In this sec- 
tion, we propose using exhaustive direct search (EDS) 
with latin hypercube sampling (LHS) as the benchmark 
method for inverse problems. It combines the accuracy 
of EDS methods, which have been widely used since the 
1960s [61], with the sampling efficiency of LHS. 

LHS [62] is a sampling technique that ensures a com- 
prehensive exploration of the solution space while main- 
taining desirable properties. The range of each variable 
was divided into equally probable intervals, and the sam- 
ple points were strategically placed to satisfy the require- 
ments of a Latin hypercube, where each sample was the 
only one in its corresponding axis-aligned hyperplane. 
The advantage of LHS lies in its independence from the 


number of dimensions, which allows for efficient sam- 
pling without the need for an increasing number of sam- 
ples. Furthermore, LHS facilitates sequential sampling, 
enabling the tracking of the already selected samples. 


B. Exhaustive Direct Search with 2 Steps Latin 
Hypercube Sampling (LHS2STEPS) 


To enhance the benchmark algorithm and mitigate the 
uncertainties introduced by the randomness of the sam- 
pling process, [13] introduced a two-step sampling strat- 
egy. This strategy involves further sampling the neigh- 
borhoods of the ту best candidates generated during the 
initial stages of the LHS algorithm. For the detailed 
pseudocode, please refer to 1, excluding the part related 
to ALHS. 


C. Exhaustive Direct Search with Assembled 2 Steps 
Latin Hypercube Sampling (ALHS) 


To improve our benchmark algorithm and further re- 
duce the uncertainty caused by the randomness of the 
sampling process, we introduce a mean mechanism at 
the end of the classical LHS algorithm to calculate the 
mean of the nı top-rated candidates in the five clones, 
which are spread around the nz top-rated candidates in 
the first sampling process. The concept is illustrated in 
Algorithm 1. 


D. Efficient Deterministic Machine Learning Algorithm: 
KNN 


In addition to the benchmark LHS method, another ef- 
ficient deterministic algorithm that can be considered for 
solving inverse models is the KNN algorithm. KNN is 
a nonparametric classification and regression algorithm 
that can be adapted to solve inverse problems. 

The KNN algorithm operates based on the principle of 
similarity. Given a new input data point, KNN finds the 
K-nearest neighbors in the training dataset based on the 
Euclidean distance and then generates the output value 
by averaging the output values of the K-nearest neigh- 
bors. The selection of the optimal value of K can be 
determined by graphing the relative Lə error between 
the predicted observation vector and the true vector, as 
shown in Figure 2(b). The analysis in Figure 2(b) reveals 
that selecting K = 1 yields the optimal outcome in terms 
of minimizing the Lə prediction error, thereby emphasiz- 
ing its advantageous accuracy. 


E. Metaheuristics Algorithms with Physical Constraints 


Recently, the application of metaheuristic algorithms 
with physical constraints [63—65] has gained significant 
attention for inverse modeling. These algorithms offer 
powerful gradient-free optimization techniques that can 
effectively handle complex inverse problems while sat- 
isfying the physical constraints imposed by system dy- 
namics. 

In this manner, we implemented and adapted practi- 
cal algorithms and constructed related solvers that were 
integrated into the developed RODT framework. These 
solvers are designed to address the inverse problem pre- 
sented by Equation 2.12. 


1. Differential Evolution Algorithm (DE) 


The differential evolution (DE) algorithm [66] offers a 
promising approach for tackling inverse modeling prob- 
lems. By leveraging its ability to handle nonlinear and 
nonconvex objective functions, DE can effectively ex- 
plore the parameter space. Furthermore, DE provides 
mechanisms to incorporate physical constraints, ensur- 
ing that the generated solutions adhere to system dynam- 
ics. The concrete procedure of the DE algorithm can 
be found in a previous work [13], in which we used the 
best mutation strategy recommended in [67] for the mu- 
tation step in DE, and the mutation donor V; с is given 
by Equation 3.14. 


Vic = Hrs,G + H« (B, G — Hra,G + Mrs,G — Hra,G), 
(3.14) 
where rz c, k € {1, 2, 3, 4,5} are randomly chosen from 
(1,2, ... NP, NP is the number of populations, и is a 
parameter, and F' and G denote the scaling factor and the 
number of generations, respectively. 


2. Particle Swarm Optimization Algorithm (PSO) 


Particle swarm optimization (PSO) [68] is a 
population-based optimization technique inspired 
by the collective behavior of bird flocking or fish school- 
ing. It involves particles representing potential solutions 
that update their positions based on the personal and 
global best positions. The algorithm iterates a set 
number of times by adjusting the particle velocities and 
positions. The best solution, represented by the global 
best position with the lowest objective function value, is 
returned. 


Algorithm 1 Exhaustive Direct Search with Assembled 2 Steps Latin Hypercube Sampling 


Inputs: 
Observations: Y, 
Initial guess: Шиа 
Forward function: ЕМ" 
Transformation operator: H 
First stage sampling number: n °! 
Second stage sampling number: n? 
U* = LHS (number = n°", initial = ръка) 
pi = argmin (|| HV ЕМЕ (u^) - Yolla) 

C]. n? 


for i from 2 to nı do 
М! = LHS ( number = n?', initial = шша) 
pi = argmin (НУ F^" (u^) - Yolla) 


end for 
for j from 1 to n; do 
и? = LHS ( number = n?? the initial = už) 
end for 
for k from 1 to n2 do 
ma argmin (|| EV, FM (и) — Yolla) 
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EUU; 
end for 
u= [Xe и) /10 
Outputs: 


Parameter: m 
Reconstructed field: V, MU (и*) 


> nı = 5 in this article. 


> Start of the main part of ALHS, n» = 10 in this article. 


> End of the main part of ALHS. 


3. Fast Simulated Annealing Algorithm (FSA) 


The simulated annealing (SA) [69] algorithm, which 
was initially proposed for solving the well-known travel- 
ling salesman problem (TSP) [70], is a popular approach 
for solving optimization problems in various fields [71], 
including nuclear engineering [72]. 

SA was inspired by metallurgical annealing, in which 
a metal is heated and gradually cooled to achieve a more 
ordered configuration. Similarly, the SA starts from 
a high-energy state (initial solution) and progressively 
lowers the temperature until it converges to a state of 
minimal energy (optimal solution). 

The fast simulated annealing (FSA) [73] is an adapta- 
tion of classical simulated annealing (CSA) [69], which 
improves the computational efficiency and convergence 
speed. CSA utilizes a local sampling approach and con- 
trols the fluctuation variance through an artificial cooling 
temperature, denoted as T.(t). FSA modifies the cool- 
ing schedule to accelerate convergence using a cooling 
temperature T(t) that decreases reciprocally rather than 
logarithmically with time t. 


4. Cuckoo Search Algorithm (CS) 


In the cuckoo search (CS) algorithm [74], nests are 
used to represent potential solutions, and eggs within the 
nests symbolize candidate solutions. Through an itera- 
tive process, the algorithm generates new solutions by 


employing Lévy flights and subsequently evaluates their 
fitness. If a newly generated solution is superior to an 
existing one within a given nest, it replaces the incum- 
bent solution. To maintain population diversity and im- 
prove the overall quality of solutions, the population is 
sorted based on fitness, and a predetermined number of 
solutions are replaced with new solutions generated via 
random walks. The algorithm terminates upon reaching 
a predefined maximum number of generations, and the 
best solution, characterized by the lowest objective func- 
tion value, is returned as the final outcome. 

CS utilizes a balanced combination of a local random 
walk and a global explorative random walk, controlled 
by a switching parameter pa. A local random walk is 
defined as follows: 


ui = pi +a H (pa — е) @ (uj — шь). (3.15) 
In the above equation, и and ju), represent two solutions 
randomly selected through permutations. H(e) denotes 
the Heaviside function, & represents point-to-point mul- 
tiplication, € is a random number drawn from a uniform 
distribution, and s represents the step size. The global 
random walk, on the other hand, is conducted using Lévy 
flights, as follows: 


ui. = ш aL(s, X), (3.16) 


where L(s, А) is defined as 


AP(A)sin(rA/2) 1 
gl tA’ 


L(s, À) (s > so > 0). 


(3.17) 


T 


5. Artificial Neural Network (NN) 


The neural network (NN) algorithm [75], inspired by 
the behavior of the human brain, simulates intercon- 
nected artificial neurons to process information, learn 
from data, and estimate unknown parameters or states. 
They consist of interconnected artificial neurons orga- 
nized into layers that receive input signals, process 
them through weighted connections, and produce out- 
put signals using activation functions. NNs can learn 
from input-output examples by adjusting the connection 
weights during training, enabling them to capture under- 
lying patterns and correlations, making them suited for 
solving physical problems [76]. 

The combination of the Adam and L-BFGS optimiz- 
ers in an alternating manner has demonstrated efficacy in 
training NNs for complex problem solutions [77]. This 
approach leverages the strengths of the two optimizers to 
enhance training and improve the overall network perfor- 
mance. 

However, owing to NN’s continuous property, it can- 
not handle discrete problems well, particularly when the 
test data are noisy, which is further discussed in Section 
V. 


F. Hybrids of Optimization Methods for Inverse 
Problem: Balancing Coarse-Fine-Grid Search 


To address the challenges posed by the discontinuity 
of KNN surrogate forward mapping and noisy observa- 
tions, we propose hybrid algorithms that combine the 
strengths of deterministic machine-learning methods for 
coarse-grid search and metaheuristic algorithms for fine- 
grid search. This methodology can also be found in fields 
that encompass feature extraction of aerial data [78], au- 
tomatic free-form optics design [79], and hyperparame- 
ter optimization [80]. The objective of a coarse-fine-grid 
search is to strike a balance between global exploration 
and local exploitation [81] to achieve accurate, efficient, 
and less sensitive parameter estimation. The hybrid op- 
timization hub designed accordingly can be further illus- 
trated by the following schema: 


1. Cuckoo Search with Differential Evolution (CSDE) 


In the later stages of the standard CS process, notable 
issues arise in the form of information waste and slow 
convergence. These issues stem from the execution of 


10 


independent evaluations and the lack of sufficient mech- 
anisms for information sharing within the population. 
Consequently, valuable information fails to propagate ef- 
fectively, resulting in suboptimal convergence rates and 
reduced search efficiency. 

From an alternative perspective, the underlying bio- 
logical nature of mutations within the standard DE algo- 
rithm was explored. In this context, a novel mathemat- 
ical elucidation of these mutations is sought. Notably, 
researchers in the natural sciences have observed that 
mutations can be understood as sequences of transfor- 
mations that transpire within a DNA molecule. Further- 
more, the transformations pertaining to modifications ex- 
hibited a statistical distribution that adhered to the Levy 
distribution. This finding suggests an analogy between 
mutations and Levy flight [82]. Drawing inspiration from 
this observation, it is conceivable to integrate the Levy 
flight process into the generation of the DE population. 

In this regard, we developed the cuckoo search dif- 
ferential evolution (CSDE) algorithm, which lies within 
the framework of the CS algorithm, integrating the DE 
mutation operation into the production of the population 
for each generation. The first step of CSDE is the ini- 
tialization: Define the objective function f(e). We set 
the population size N, problem dimension D, discovery 
probability Pa, step size a, scaling factor F’, crossover 
probability CR, maximum number of iterations T', and 
search domain range [juy, Hub]. The steps are presented 
in Algorithm 2. 


2. Advanced CSDE (ACSDE) and Advanced DE (ADE) 


We utilize the SVC algorithm to classify the test pa- 
rameters and denote the center of the classified hyper- 
cube as the initial machine learning estimate sv c; then, 
we use the CSDE and DE algorithms to further explore 
the parameter space with the initial population spread 
around the center usy c. In the SVC process, according 
to a previous research [15] on the importance rate of the 
input parameters, the last parameter, temperature Тп, 
contributes slightly, and we also consider that Ттт varies 
at a relatively small boundary; thus, we only administer 
the SVC process to the first three parameters. We chose 
two representative metrics to evaluate the performance of 
the SVC model. 

Confusion Matrix. The confusion matrix [83] pro- 
vides a comprehensive overview of the binary classifi- 
cation results by displaying the number of true positive 
(TP), true negative (TN), false positive (FP), and false 
negative (FN) predictions. This allowed us to assess 
the accuracy of the model’ by differentiating classes and 
identifying potential sources of misclassification. We can 
infer from Figure 5 that our SVC classification model 
predicts the first two parameters with relatively high ac- 
curacy while exhibiting confusion concerning the third 


Coarse-grid search Fine-grid search 


KNN 


LHS 


SVC 


LHS 


DE 


CSDE 


Hybrid optimization hub 


Fig. 4: Hybrid optimization schema of coarse-fine-grid search. 


parameter. 

Fl-score. The Fl-score [84], defined in Equation 
3.18, is a popular performance metric for binary clas- 
sifier models. In Equation 3.18, ТР, FN, FP are the 
numbers of true positives, false negatives, and false pos- 
itives classified by the model, respectively. 


Е ТР 
ТР+і(ЕР + FN) 


The predictive model employed in this study for each 
parameter is a multi-class SVC utilizing the “one-vs- 
others” strategy. This approach involves training sep- 
arate binary classifiers to distinguish between the can- 
didate intervals for a given parameter. Consequently, 
performance evaluation necessitates the adoption of the 
Macro Fl-score as a metric. The Macro Fl-score 3.19 
is computed as the average of the Fl-scores obtained for 
each interval, reflecting the overall effectiveness of the 
classifiers in capturing intra-parameter variations: 


1 N 
Flinacro :— N 3 F1;, 


where 7 is the interval index, and № is the number of 
intervals. 

Based on the Confusion Matrix shown in Figure 5 and 
F lmacro, We designed the grid size of the fine-grid search 
process as follows: 


Fi (3.18) 


(3.19) 


3. K-Nearest Neighbour plus Exhaustive Direct Search with 
Latin Hypercube Sampling (KNNLHS) 


Integrating KNN with LHS provides a robust method- 
ology to address inverse problems. To enhance the com- 
prehensiveness of our research, we standardized the grid 
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size for a coarse-fine grid search. This standardization is 
consistent with the specifications listed in Table 1. 

In conclusion, the k-nearest neighbour plus exhaustive 
direct search with latin hypercube sampling (KNNLHS) 
method is a robust and efficient approach compared 
to similar strategies. This technique combines the ex- 
ploratory capabilities of KNN in conducting coarse-grid 
searches with the efficiency and diversity of LHS in sam- 
pling local solution spaces. Consequently, the KNNLHS 
facilitates precise and efficient parameter estimation, ac- 
curate navigation of complex objective functions, and the 
management of noisy datasets. The effectiveness of this 
approach and its comparative advantages are discussed 
in Section V. 


IV. APPLICATION TO NUCLEAR REACTOR 
A. Setup of Reactor Physics Operation Problems 


Setup of background information for HPR1000. 
The objective of constructing the RODT was to accu- 
rately predict the power distribution within the HPR1000 
reactor core during operation. The core consists of 
177 vertical nuclear fuel assemblies, 44 of which are 
equipped with self-powered neutron detectors (SPNDs) 
to measure neutronic activity and power fields. Figure 
6 provides a visualization of a horizontal slice of the 
HPR1000 core and an axial slice of an SPND-equipped 
assembly, with only one-quarter of the core displayed 
owing to the symmetry along the z and y axes. The gray 
fuel assemblies represent those containing control rods, 
whereas the assemblies marked with D indicate the pres- 
ence of SPNDs. For more detailed information regard- 
ing the HPR1000 reactor and generic neutronic physi- 


Algorithm 2 Cuckoo Search Algorithm with Differential Evolution 


Step 1: Population initialization 
Initialize the population Р and the initial position using a random distribution within the search domain range [ui Hub] 
while С < Gmax do 
for each individual и; in population Рс do 

Step 2: Population updating 

Generate random indices r1, r2, r3 (different from i) 

for each dimension j in problem dimension dim(u) do 

if random(0, 1) « CR then 
u$; = u£; + ®® (ифз — иез) @ L (s, A) 


else 
и) = ие; 
end if 
end for 
Step 3 Mutation Step 


Generate random indices r3, r4, r5 (different from i) 
Generate a donor vector V^ € Rm: 
М = us БЕ (um — Ша) 
Step 4 Crossover Step 
Generate a trial vector US є R?/mG2. 
for j = 1 to dim() do 
if random(0, 1) < Cr then 


G G 
U, ig — V; j 
else 
G G 
Ui; j = №, 
end if 
end for 
Step 5 Selection Step 


Evaluate the trial vector U, E 
if f(U7) < f(u?) then 


up = UË 
else 
Hi = п? 
end if 
Step 6 Abandon nests 


Sort the population Рс based on the fitness values in ascending order. 
Determine the number of solutions to be replaced, терасе = rand(pa · N P). 
for i = 1 to replace do 
Generate a new solution X in the search space 
with random walk given in Equation 3.15. 
Evaluate the objective function f for X. 
if f(X) < f(worst solution їп Р) then 
Replace the worst solution in Pg with X. 
end if 
end for 
end for 
G=G+1 
end while 


TABLE 1: Grid size design for the Coarse—Fine-Grid search based on the SVC prediction evaluation. 


ParameterNames St Bu Pu 

Васо 9.98E-01 9.93E-01 5.18E-01 
Size coarse-grid(SVC) 41 steps 250MWd/tU 8% FP 

Size fine-grid(DE/CSDE) 41 steps 250 MWd/tU 20% FP 
Relative Size coarse-grid(SVC) 6.67% 10.00% 10.00% 
Relative Size fine-grid(DE/CSDE) 6.67% 10.00% 25.00% 
Parameter Boundary [0,615] [0,2500] [2096, 10096] 
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((a)) 5 


((b)) Bu 


((с)) Pw 


Fig. 5: Confusion Matrix Plot that evaluates the efficacy 
of SVC prediction in the first 3 parameters (S+, Bu, Pw, 
respectively). 


cal model, please refer to [14]. Further descriptions of 
the model for data assimilation and the initial work on 
RODT can be found in [35] and [13], respectively. 
During the normal operation of the HPR1000 reactor, 
two types of control rods were utilized for regulation. 
The first type, known as compensation rods, is responsi- 
ble for coarse control and a substantial reactivity reduc- 
tion. These compensation rods comprise four subtypes: 
G1, G2, N1, and №2. The second type, called regulat- 
ing rods (R rods), are employed for fine adjustments to 
maintain the desired power or temperature [85]. Power 
evolution within the reactor is influenced by various fac- 
tors, including the movement of the control rods, the bur- 
nup of nuclear fuel, variations in the power level of the 
reactor core, and fluctuations in the inlet coolant temper- 
ature. This evolution was mathematically modeled us- 
ing two-group diffusion equations [13] and numerically 
solved using the CORCA-3D code package [86]. Devel- 
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oped at NPIC, the CORCA-3D code is capable of solving 
3D few-group diffusion equations, considering thermal- 
hydraulic feedback, and performing pin-by-pin power re- 
construction. 

Setup of parameter domain. In this study, we con- 
sider CORCA-3D to be an opaque solver in which we 
input the parameter vector yz and obtain the output Ф. In 
particular, the power field ® is constrained to be depen- 
dent on a set of "general" parameters that indicate the 
stage of the reactor's life cycle: 


H = (Si, Bu, Pu, Tin). (4.20) 


* S,: The control rod insertion steps, ranging from 
0 to 615, represent the movement of the compen- 
sation rods from all rod clusters being out (ARO) 
to fully inserted. This range takes into account the 
overlap steps. 


B,,: The average burnup of the fuel in the entire 
core indicates the amount of energy extracted from 
the fuel and increases over time. Its value ranges 
between 0 (at the beginning of the fuel’s life cycle) 
and By, max, Which is set to 2500 MWd/tU (the end 
of the fuel’s life cycle). The specific evolution of 
В depends on the reactor's operational history. 


Pw: The power level of the reactor core ranges 
from 0.3 to 1 FP (full power). 


Tin: The temperature of the core coolant at the 
inlet falls within the range of 290 to 300°C. 


In this manner, can be implicitly represented 
by Ф(и) Фф (St, Bu, Pw, Tin), thanks to CORCA- 
3D. There are 177 fuel assemblies in HPR1000, and 
each assembly is numerically represented using 28 
vertical levels. Thus, ® is a vector with dimensions 
N = 4956(— 177 x 28). The discrete solution set M = 
(é(u) € RN | we D} consists of P = 18480 solu- 
tion snapshots with the parameter configuration D : 


B; 59 PE ӘТ, where Sf = {0,1,...,615}, B5 = 
(0, 50, 100, 150, 200, 500, 1000, 1500, 
2000, 2500}, P =  RU3(30,100) and Ti = 


RU3(290,300). Тһе operator RU3(a,b) represents 
three independent and identically uniformly distributed 
samplings in the closed set [a,b] where а < b. In this 
study, we selected 9096 of 18480 snapshots for training 
the forward and inverse models, and the remaining 10% 
of snapshots were used for test purposes. 

Setup of observation data and modeling the obser- 
vation noise. In the context of RODT, the observations 
denoted by Y, are used to infer the input parameter vector 
ш and subsequently determine the resulting power field 
Ф. The arrangement of these observations, organized 
node by node, is shown in Figure 6. Each observation 
value represents the average value over the correspond- 
ing node in which the SPNDs are located. It is impor- 
tant to clarify that the observations used in the analysis 
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Fig. 6: Quarter of the core in the radial direction (white square: fuel assembly with SPNDs, grey square: fuel 
assembly with control rods, D: fuel assembly with neutron detectors) [13]. 


were not obtained from direct engineering measurements 
but were derived from numerical simulations of ® using 
CORCA-3D. 

The noise associated with the observations can be 
modeled as either a Gaussian or a uniform distribution, 
depending on the specific requirements. Specifically, the 
noise term can be mathematically represented as follows: 


e = Y, - N(0,0,Y,) or Yo - U(—0, ô, Yo) 


In this context, the vector N (0, с, Yo) denotes a collec- 
tion of random variables following a normal distribution. 
These random variables have the same dimensions as Y, 
and are characterized by a mean of zero and standard 
deviation of o. Similarly, the vector U (—ó, ó, Yo) rep- 
resents a set of random variables uniformly distributed 
within the interval [—ó, ó]. These variables have the same 
dimensions as Y,. 

Considering the prevalence of Gaussian noise in real- 
world nuclear reactors, our subsequent analysis primarily 
focused on evaluating the performance of the algorithm 
under Gaussian noise conditions. 

Software implementation We integrated our software 
into the RODT platform using the following modules: 


* rodtPro: This module implements data cleaning 
and preprocessing techniques based on domain ex- 
pertise. 


* rodtROM: This module constructs ROM through 
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techniques like POD, SVD, and RB. 


* rodtML: This module focuses on surrogate mod- 
eling and inverse model optimization. It also in- 
cludes hyperparameter optimization techniques. 


rodtPost: This is a post-evaluation module that de- 
fines convergence metrics and conducts compara- 
tive analyses using evaluation methodologies de- 
rived from peer-reviewed literature. 


rodtUI: This module is responsible for visualiz- 
ing the results, providing an intrinsic review of the 
methods, showcasing the optimized hyperparame- 
ters, and designing user-specific visualizations. 


B. Practice Metrics for Inverse Model Evaluation 


To evaluate the performance of the algorithms, we 
designed a series of numerical experiments using nu- 
clear reactor data produced by an intrusive program with 
added noise. We carefully selected metrics related to ac- 
curacy and finally defined the normalized reconstruction 
fields Lə error (Es) and La error (Eso) to evaluate the 
performance of the comparative algorithms. 

The reconstructed field Lz error is defined in Equation 
4.2] to evaluate the average estimation of the field, where 


1/2 
ll e ||2 := b» (e?) denotes the L;norm. 


E» = |Q FM (и*) = Фе |l / Il Фе lla 


The reconstructed field Lo error is defined in Equa- 
tion 4.22 to evaluate the worst field estimation, where 
|| e ||. := max; |e;| denotes the Lo; norm. 


Ev := ||QF™™ (и*) — Ф. || / Фе llo (422 


It should also be noted that the reconstructed observa- 
tion Lə error and La can be generated similarly. 


(4.21) 


V. NUMERICAL RESULTS 


In this section, we describe the experiments conducted 
to evaluate the performance of parameter identification 
and state estimation for the proposed algorithms. The 
evaluation was performed on a test dataset consisting of 
40 specimens from HPR1000, as described in Section IV. 
All the proposed inverse modeling approaches are uti- 
lized in the experiments, and it is important to mention 
that the cost function employed remains consistent with 
Equation 2.12. Notably, we focused solely on testing the 
CSDE method and did not include the CS method. Fur- 
thermore, we introduce the proposed accuracy metric, il- 
lustrated in Section IV B, to assess its performance. 


A. Accuracy and Stability Analysis of the Parameter 
Identification Phase 


In this section, we describe the tests conducted to eval- 
uate the accuracy of parameter identification. Initially, 
we assessed the effectiveness of the proposed hybrid 
methodology, the coarse-fine-grid search, as discussed 
in Section ШЕ, using the KNNLHS method as an illus- 
trative example. To visualize this process, we generated 
scatter plots depicting the predicted values of the param- 
eters 5; and P,,, which were identified as the two dom- 
inant parameters in a previous study [15] (see Figure 7 
for six cases with different input parameter settings. 

Figure 7 shows a visualization of the local search 
space in the inverse algorithms. The figures in the first 
row depict scenarios with varying P, values ranging 
from relatively low to moderate to relatively high. Corre- 
spondingly, the figures in the second row illustrate cases 
with varying 5; values, again encompassing relatively 
low, moderate, and high values. These figures effectively 
show the considerable improvements achieved using the 
KNNLHS method, particularly in the simultaneous op- 
timization of both P,, and S;. The visualizations un- 
derscore the effectiveness of the method in yielding en- 
hanced optimization outcomes. 

The parameter prediction deviation is demonstrated in 
Figures 8 over the test set. Owing to space limitations, 
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we chose only KNN, KNNLHS, and FSA for illustra- 
tion and tested clean and noisy (с = 5%) observation 
data. Using the LHS algorithm, the KNN prediction for 
parameter identification became more accurate, and the 
prediction deviation for the FSA algorithm was concen- 
trated at relatively low intervals, thus proving its excel- 
lent accuracy. In general, the accuracy of burnup, the 
power level, and the control-rod step are acceptable from 
an engineering perspective. 


B. Comparison of the Convergence Rate among 
Metaheuristic Algorithms 


The convergence rate is a crucial criterion for eval- 
uating algorithms that are integrated with the genera- 
tion design; thus, in this section, we compare different 
metaheuristic algorithms and their hybrids. Their con- 
vergence performance in 50 iterations (generation) with 
clean and noisy (c = 1% and с = 5%) observation data 
is shown in Figure 9. 

When the observation data were clean, the FSA ex- 
hibited the fastest convergence rate, followed by the ad- 
vanced differential evolution (ADE). Moreover, when in- 
tegrated within an SVC coarse-grid search, both the DE 
and CSDE algorithms demonstrated the ability to escape 
the local optima in the early stages while incurring a rel- 
atively low-cost function, thereby surpassing their stan- 
dalone performance. 

Among the two population-based algorithms, DE out- 
performed CSDE. This can be attributed to two fac- 
tors: (i) the adoption of the best mutation method for 
DE, where we use the optimal mutation strategy for DE 
while using the classical strategy; and (11) the relatively 
straightforward nature of our problem, wherein the in- 
verse problem lacks multiple optima. Consequently, the 
CS step in each generation of CSDE may lead to unnec- 
essary exploration of the parameter space that is distant 
from the optimal solution. 

In the presence of noise-contaminated observation 
data, FSA continues to exhibit superior convergence ca- 
pabilities. However, its advantage diminishes as the level 
of Gaussian noise reaches 5%, eventually leading to a 
decline in advantage. Notably, the ADE and advanced 
cuckoo search differential evolution (ACSDE) failed to 
deliver satisfactory results under noisy conditions. This 
is primarily due to the erroneous predictions of SVC, 
which further leads to inadequate initial estimates pro- 
vided by population-based algorithms. 


C. Evaluation of State Estimation Phase by Accuracy 


In this section, we evaluate various methods in the 
state-estimation phase based on accuracy. First, we com- 
pared the accuracy and stability under two metrics: the 
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Fig. 7: Distribution of generated candidate parameters и! = S; and и? = P,, with hybrid algorithms. The solid 
green points are the true values, while the hollow yellow points stand for the candidate parameter generated by the 
mere KNN algorithm. The hollow blue triangles represent the five optimal samplings by LHS around the initial guess 
given by the KNN evaluated in the observation space. The solid red points represent the output of the algorithms, 
which is the best among the five optimal samplings. 


relative error in the Lz norm defined by Equation 4.21 
and the Loo norm defined by Equation 4.22 for field ®, 
which is presented in Table 2. The ‘Mean’ values in 
the tables represent the average accuracy across 40 tested 
specimens, while the ‘STD’ values indicate the standard 
deviation of errors over the same set of specimens. The 
optimal values are indicated in bold. 

In general, integration of the LHS algorithm enhanced 
the accuracy of the KNN model. This improvement can 
be attributed to the LHS sampling approach, which mit- 
igates the discontinuity issue associated with KNN pre- 
dictions by sampling around the initial estimate. 

Furthermore, the NN method exhibited a decline in ac- 
curacy when confronted with noise-contaminated obser- 
vational backgrounds. This is attributed to its reliance on 
gradients, rendering it less effective in discrete cases. 

To visualize the performance of the inverse algorithms 
directly, we randomly selected one specimen from a test 
dataset. Figures (11 and 12) depict the relative Lo error 
of the field prediction at the 11th vertical level. Readers 
are referred to Figure 10(a) for the true field values. 

To fully demonstrate the reconstruction accuracy of 
our approach, we examined cases in which the recon- 
struction was relatively poor, particularly around the con- 
trol rod interfaces. Reconstruction of the physical fields 
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in locations near the control rod interfaces is challeng- 
ing because the behavior in these regions is not easily 
captured with high fidelity. Readers can refer to Figure 
10(b) for the true field value on this axis, and Figures (13 
and 14) depict the prediction errors of various methods 
on this axis. 

However, our results also show that even in such dif- 
ficult cases, our method can achieve good reconstruc- 
tion in most cases. Although some techniques performed 
worse around the interfaces, the majority of our proposed 
methods still facilitated a reconstruction quality compa- 
rable to other areas across the domain. Overall, this 
validation indicates that the proposed approaches hold 
promise for representing complex reactor physics behav- 
iors with suitable accuracy, even when challenges arise 
locally because of design complexities, such as control 
rod insertion points. 


D. Evaluation of State Estimation Phase by Time Cost 


In this section, we present an overview of the time 
costs (we take the mean of the 40 specimens) associ- 
ated with each proposed algorithm in the state-estimation 
phase, which can be found in Table 3. The optimal val- 
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TABLE 2: Comparison of reconstructed power field errors using different algorithms. 


noise level 0 0.01 0.02 0.03 0.04 0.05 
Mean L5 

KNN 1.0662E-02 1.1128E-02 1.1461E-02 1.1461E-02 1.3296E-02 1.3541E-02 
LHS 4.0457E-02 4.0210E-02 4.0823E-02 4.0754E-02 4.1138E-02 4.1790E-02 
ALHS 3.4294E-02 3.4510E-02 3.4837E-02 3.4865E-02 3.5218E-02 3.7221E-02 
DE 5.3219E-03 6.6446E-03 7.8485E-03 9.2407E-03 1.1051E-02 1.3777E-02 
ADE 3.8119E-03 6.8077E-03 7.8185E-03 9.1585E-03 1.1134E-02 1.3526E-02 
CSDE 1.7029E-02 1.6581E-02 1.7011E-02 1.7390E-02 1.8874E-02 2.0011E-02 
ACSDE 7.6471Е-03 4.4970E-02 4.5157E-02 4.6044E-02 4.7066E-02 4.8563E-02 
PSO 1.0093E-02 1.1502E-02 1.3355E-02 1.4640E-02 1.6516E-02 1.7939E-02 
FSA 3.8796E-03 5.1561E-03 6.5906E-03 8.1864E-03 1.0106E-02 1.2991E-02 
NN 2.6737E-02 9.8060E-02 1.3113E-01 1.6029E-01 1.8785E-01 2.1606E-01 
LHS2STEPS 3.1108E-02 3.1410E-02 3.1664E-02 3.2919E-02 3.3712E-02 3.4975E-02 
KNNLHS  7.0615E-03 9.0129E-03 9.6693E-03 1.0514E-02 1.2406E-02 1.3408E-02 


STD L2 
KNN 5.3359E-03 5.5301E-03 5.2201E-03 5.2201E-03 5.5587E-03 5.7721E-03 
LHS 2.5494E-02 2.5404E-02 2.5938E-02 2.5916E-02 2.5778E-02 2.5591E-02 
ALHS 2.3719E-02 2.3621E-02 2.3377E-02 2.3131E-02 2.2688E-02 2.2074E-02 
DE 5.3195E-03 4.6991E-03 4.6245E-03 4.6625E-03 5.7638E-03 8.1052E-03 
ADE 3.1455E-03 4.6463E-03 4.5439E-03 4.6656E-03 6.0923E-03 7.7264E-03 
CSDE 9.7295E-03 7.3634E-03 7.3791E-03 7.2028E-03 7.1571E-03 7.0336E-03 
ACSDE 5.2449Е-03 2.2811E-02 2.4206E-02 2.4100E-02 2.3438E-02 2.2736E-02 
PSO 1.2374E-02 5.3048E-03 6.9606E-03 6.2216E-03 6.6237E-03 7.2507E-03 
FSA 6.2602E-03 3.5832E-03 3.5640E-03 3.7438E-03 4.6444E-03 6.0076E-03 
NN 3.4667E-02 7.3152E-02 8.3294E-02 9.8081E-02 1.1251E-01 1.2140E-01 
LHS2STEPS 2.1742E-02 2.1617E-02 2.1475E-02 2.1528E-02 2.1408E-02 2.0572E-02 
KNNLHS  4.9637E-03 5.3025E-03 4.9676E-03 5.8222E-03 5.9288E-03 6.1111E-03 


Mean La; 
KNN 4.1213E-02 4.3991E-02 4.4840E-02 4.4840E-02 4.9788E-02 4,8775E-02 
LHS 1.3126Е-01 1.2567E-01 1.2326E-01 1.1998E-01 1.1904E-01 1,2479E-01 
ALHS 1.0932E-01 1.1062E-01 1.1454E-01 1.1529E-01 1.1662E-01 1,2263E-01 
DE 2.4302E-02 2.7539E-02 3.2492E-02 3.8239E-02 4.4672E-02 5,2324E-02 
ADE 1.5357E-02 2.8079E-02 3.2408E-02 3.7967E-02 4.4269E-02 5,1787E-02 
CSDE 5.8757E-02 6.0849E-02 6.3653E-02 6.5595E-02 6.9252E-02 7,2312E-02 
ACSDE 3.1167Е-02 1.4168E-01 1.4022E-01 1.4177Е-01 1.4747E-01 1,4928E-01 
PSO 4.6008E-02 4.6052E-02 5.2500E-02 5.6560E-02 6.2328E-02 6,7784E-02 
FSA 1.6543E-02 2.2407E-02 2.7837E-02 3.3938E-02 4.0081E-02 4,9525E-02 
NN 8.1834E-02 2.3886E-01 3.0443E-01 3.6030E-01 3.8488E-01 4,1158E-01 
LHS2STEPS 1.0144E-01 1.0014E-01 1.0167E-01 1.0810E-01 1.1415E-01 1,1607E-01 
KNNLHS 2.7072Е-02 3.7594E-02 4.0961E-02 4.3185E-02 4.7107E-02 5,0309E-02 


STD Loo 
KNN 2.7190E-02 3.0633E-02 2.8849E-02 2.8849E-02 2.8941E-02 2.9448E-02 
LHS 8.7321E-02 8.2894E-02 7.8177E-02 7.6610E-02 7.7353E-02 8.1824E-02 
ALHS 8.8403E-02 8.9334E-02 8.8679E-02 8.8367E-02 8.9894E-02 8.7995E-02 
DE 2.2486E-02 1.8169E-02 1.9161E-02 1.9406E-02 2.4870E-02 3.2087E-02 
ADE 1.5143E-02 1.8220E-02 1.9221E-02 2.0268E-02 2.5709E-02 3.0895E-02 
CSDE 4.0726E-02 2.6558E-02 2.8031E-02 2.8521E-02 2.6937E-02 2.9098E-02 
ACSDE 2.1241Е-02 6.5362E-02 6.7719E-02 7.0399E-02 6.9278E-02 6.8662Е-02 
PSO 5.6853E-02 2.4197E-02 2.5349E-02 2.5845E-02 2.8182E-02 3.1462E-02 
FSA 2.3181E-02 1.5810E-02 1.5634E-02 1.6278E-02 2.0246E-02 2.6014E-02 
NN 9.1032E-02 1.6232E-01 1.8105E-01 1.8648E-01 2.0080E-01 1.9789E-01 
LHS2STEPS 7.1699E-02 6.9753E-02 7.1235E-02 7.3554E-02 7.4717E-02 7.2363E-02 
KNNLHS  2.2045E-02 2.6180E-02 2.4270E-02 2.8580E-02 2.9683E-02 3.2208E-02 


18 


Iteration 


((a)) Clean Observation 


Iteration 


(b) с = 1% 


3x107} 


2x107} 


6x107? 


Iteration 


((c)) c = 5% 


Fig. 9: Convergence rate of the cost function (defined by 
the Lo norm) with different metaheuristic optimizers. 
The x-axis represents the iteration (generation), and we 
set the maximum iteration (generation) to 50. 


ues are highlighted in bold, and the second optimal value 
is underlined. 

It is worth mentioning that the NN method demon- 
strated notable speed. However, metaheuristic algo- 
rithms incur relatively high time costs. This is primarily 
attributed to the repeated computation of the cost func- 
tion in Equation 2.12, which occurs 50 times for each 
population member over 50 iterations/generation cycles, 
resulting in a total of 2500 cost function evaluations. 
Moreover, these algorithms involve additional compu- 
tations and sorting operations for information exchange 
and identification of the best solutions, contributing to 
the discrepancy in time cost. 

Furthermore, it is important to note that the meth- 


19 


1 

2 
3 14 

4 
5 12 

6 
1.0 

7 
8 0.8 

9 
0.6 

10 
n 0.4 

12 
43 0.2 
14 0.0 

15 


ABCDEFGH tI J KLMNO 


((a)) Radial power distribution of the 11th axial plane 


16 


10 15 25 
Axis Z 


20 


((b)) Axial power distribution of ће DO1 fuel assembly 


Fig. 10: True value of the power distribution on one 
specimen in the test dataset over the core in a realistic 
HPR1000 reactor. 


ods LHS, ALHS, and LHS2STEPS have approximately 
the same time costs. In the algorithm LHS, we set the 
number of samples to 1000. In algorithms ALHS and 
LHS2STEPS, we set the first-stage sampling number 
to 750, and the second-stage sampling is repeated five 
times, with each repetition sampling number uniformly 
set to 50, adding up to 1000 in total. Model-based al- 
gorithms, such as KNN and its hybrid KNNLHS, have 
relatively low time costs, that is, less than one second, 
which is acceptable in the industry. 


E. Comprehensive Evaluation of Various Algorithms 


In the state-estimation phase, a comparison of the pro- 
posed algorithms in terms of time cost and accuracy is 
presented in Figure 15. The KNNLHS algorithm demon- 
strates a commendable level of accuracy, achieving solu- 
tions within 0.1 s for both clean and noise-contaminated 
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Fig. 11: Reconstructed relative errors of radial power distribution on the 11th axial plane with clean observations. 


observation backgrounds. Consequently, KNNLHS can solutions. 
be considered suitable for real-time online inverse mod- 
eling in the context of RODT. Moreover, in scenarios in 
which time constraints are not limiting factors, the FSA, 
DE, and ADE algorithms provide reasonably accurate 
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Fig. 12: Reconstructed relative errors of radial power distribution on the 11th axial plane with observations of noise 
level с = 5%. 


VI. CONCLUSION AND FUTURE WORKS 


In this study, we developed and implemented a series 
of powerful algorithms for the RODT, aiming to address 
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the challenges in parameter identification and state esti- 
mation for nuclear reactor monitoring and optimization 
applications. Novel hybrid optimization algorithms, in- 
cluding ADE, CSDE, ACSDE, and KNNLHS, were de- 
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Fig. 13: Predicted relative error for the DO1 fuel assembly with clean observations. 


TABLE 3: Comparison of time cost for different algorithms. 


noise level 0 0.01 0.02 0.03 0.04 0.05 
KNN 6.3696E-03 3.1396E-03 3.1599E-03 3.1594E-03 3.1915E-03 3.1527E-03 
LHS 7.3946E-01 7.3326E-01 7.2510E-01 7.3279E-01 7.2927E-01 7.2180E-01 
ALHS 7.1437E-01 7.3674E-01 7.1443E-01 7.1613E-01 7.1730E-01 7.1050E-01 

DE 3.0745Е+00 2.9459E+00 2.8510E+00 2.9223E+00 2.8299E+00 2.7874E+00 
ADE 2.9249E+00 2.9339E+00 2.8915E+00 2.9284E+00 2.9743E+00 3.0290E+00 
CSDE 4.7069E+00 4.6679E+00 4.6691E+00 4.6531E+00 4.6615E+00 4.6227E+00 

ACSDE = 4.8054E+00 4.6563E+00 4.6374E+00 4.6886E+00 4.6303E+00 4.7069E+00 

PSO 2.4404E+00 2.5169E+00 2.4969E+00 2.3198E+00 2.3551E+00 2.3151E+00 
FSA 2.7848E+00 2.8079E+00 2.7324E+00 2.7212E+00 2.7090E+00 2.6977E+00 

NN 4.1338E-04 4.2277E-04 4.0678E-04 4.0900E-04 4.0602E-04 4.0516E-04 

LHS2STEPS 7.0137E-01 7.0546E-01 7.0171E-01 7.0508E-01 7.0827E-01 7.0630E-01 

KNNLHS  7.0515E-02 6.4439E-02 6.3790E-02 6.4109E-02 6.3136E-02 6.3097E-02 
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Fig. 14: Predicted relative error for the D01 fuel assembly with observations of noise level с = 5%. 


signed and implemented to effectively solve the inverse 
problems posed by the discontinuous surrogate forward 
modeling approach. Both deterministic and metaheuris- 
tic optimization methods have been investigated for their 
applicability to nuclear engineering systems character- 
ized by high-dimensional, nonlinear relationships. The 
integration of coarse and fine grid searching strikes an 
optimal balance between computational efficiency and 
solution accuracy. 

Through extensive experimentation and evaluation, 
we conducted a thorough comparative analysis of param- 
eter identification accuracy, field reconstruction accu- 
racy, and convergence profile. FSA exhibited the fastest 


23 


convergence rate, and the proposed ADE method fol- 
lows closely behind. Moreover, integrating standard DE 
or CSDE algorithms within the SVC coarse grid search 
framework enables these techniques to effectively escape 
local optima in early iterations while maintaining rel- 
atively low-cost function values, surpassing their stan- 
dalone implementations. Additionally, it is observed 
that pairing the LHS approach within the KNN surro- 
gate model enhances the accuracy by mitigating discon- 
tinuity issues through local sampling around the initial 
guess. Thus, our approach effectively addresses the chal- 
lenges of discontinuity and non-differentiability in sur- 
rogate forward mapping generated by the KNN model. 
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Fig. 15: Comparison of the synthetic performance for 
different inverse algorithms. 


Conversely, the NN method experiences declining ac- 
curacy when faced with noise-contaminated measure- 
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ments, likely owing to overreliance on gradient informa- 
tion, which proves less effective under such uncertainty. 

Our algorithms demonstrated high computational effi- 
ciency, with solutions obtained within | s and as fast as 
< 0.08 s for KNN, KNNLHS, and NN. They also ex- 
hibited robustness in handling noise-contaminated back- 
grounds. Our work contributes to the realization of 
RODT by providing effective, efficient, and reliable sup- 
port for real-world nuclear power plant operation and 
sensor data. These algorithms establish a promising 
foundation for advanced system modeling and engineer- 
ing decision support. 

By specifically tailoring solutions to nuclear inverse 
and surrogate modeling challenges, this research ad- 
vances the state-of-the-art RODT methodology. The im- 
proved parameter identification and state estimation fa- 
cilitate the optimization of reactor monitoring, control, 
and performance through data-driven DT applications. 
In the future, novel methods can contribute to our RODT 
system, and with the same methodology, we can sup- 
port the development and implementation of these trans- 
formative technologies in other industries requiring real- 
time DTs. 

Continued work validating diverse systems and pro- 
gressively integrating multi-physics modeling holds 
promise for fully realizing the potential of RODTs in 
transforming nuclear engineering practice with power- 
ful real-time and decision-support capabilities. Realizing 
this potential requires ongoing progress, collaboration, 
and open innovation across expert areas. 
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